Digital Filtering - Using filtfilt
Difficulty Level:
Tags pre-process☁filter☁filtfilt

As it has been mentioned in a previous Jupyter Notebook concerning digital filtering , recorded signals usually contain noise , which may have different origins, that has to be minimised as much as possible to obtain a high quality signal with the highest achievable signal to noise ratio.

This Jupyter Notebook will focus on the use of a specific type of digital filter: the filtfilt , also known as zero phase filtering or forward and backwards digital filter. It applies the same filter twice: once forwards and once backwards. This combination results in a filter with zero phase distortion and an order that is twice that of the original filter.

This type of filtering is mostly useful to preserve features in a filtered time waveform exactly where they occur in the unfiltered signal.

1 - Importation of the required packages

In Python, it is very easy to use public code to help our signal processing tasks. In this case, besides the biosignalsnotebooks package, we will also use some functions of the numpy Python package.

In [2]:
# biosignalsnotebooks own package for loading and plotting the acquired data
import biosignalsnotebooks as bsnb

# Scientific packages
from numpy import arange, sin, pi
from numpy.random import randn
In [2]:
# Base packages used in OpenSignals Tools Notebooks for plotting data
from bokeh.plotting import figure, output_file, show, curdoc
from bokeh.io import output_notebook
from bokeh.layouts import column, row
from bokeh.models import ColumnDataSource, Plot, LinearAxis, BoxAnnotation, Arrow, VeeHead, LinearAxis, Range1d
output_notebook(hide_banner=True)

2 - Use of filtfilt in a generated sine function

In this section, we will start our analysis in a synthetic signal in order to ease its analysis. Thus, we will construct a time axis and associate a sine wave to it. Then, we will add random generated noise to the sine wave. Finally, we will apply a normal filter and a filtfilt to compare the results.

First, let"s build the sine wave:

In [3]:
# First, we will construct the time axis
sine_time = arange(0, 2, .01)
# Now, build the sine wave associated with the time axis. 
sine_signal = sin(pi * sine_time)

Then, we will generate random noise and add it to the sine wave:

In [4]:
# Let's add some random noise to our new signal
# First, generate the random noise
noise = randn(len(sine_time))*0.1
# Then, add the noise to the sine wave
noisy_signal = sine_signal + noise

Finally, we will apply two types of low-pass filter: one normal and the other using the filtfilt.

In [5]:
# Applying the normal filter with low cut frequency of 20 Hz
lfilter_signal = bsnb.lowpass(noisy_signal, 20, order=3)

# Applying the filtfilt filter with low cut frequency of 20 Hz
filfilt_signal = bsnb.lowpass(noisy_signal, 20, order=3, use_filtfilt=True)

The next plots show the results of the application of the different filters. The first shows the noisy sine wave without the application of filters. The middle figure shows the original signal with noise and the signal resulting from the application of a normal low-pass filter. The last figure shows the same, but the signal in dashed line is the result of the application of the filter using filtfilt .

In [6]:
# Function to help us plot the examples
def print_filtfilt_plots(time, signal, lfilter_signal, filtfilt_signal, legend="Original Signal"):
    p = figure(plot_width=400, plot_height=200)
    p.ygrid.grid_line_alpha=0.5
    p.line(time, signal, line_width=1, legend_label=legend)
    p = bsnb.applyOpenSignalsStyle(p)
    p.yaxis.axis_label = 'Value'
    p.xaxis.axis_label = 'time'

    q = figure(plot_width=400, plot_height=200)
    q.ygrid.grid_line_alpha=0.5
    q.line(time, signal, line_width=2, line_color='blue', legend_label=legend, line_alpha=0.5)
    q.line(time, lfilter_signal, line_width=2, line_dash='dashed', line_color='green', legend_label='Filter')
    q = bsnb.applyOpenSignalsStyle(q)
    q.yaxis.axis_label = 'Value'
    q.xaxis.axis_label = 'time'

    r = figure(plot_width=400, plot_height=200)
    r.ygrid.grid_line_alpha=0.5
    r.line(time, signal, line_width=2, line_color='blue', legend_label=legend, line_alpha=0.5)
    r.line(time, filtfilt_signal, line_width=2, line_dash='dashed', line_color='red', legend_label='Filfilt')
    r = bsnb.applyOpenSignalsStyle(r)
    r.yaxis.axis_label = 'Value'
    r.xaxis.axis_label = 'time'
    
    return p,q,r

p, q, r = print_filtfilt_plots(sine_time, noisy_signal, lfilter_signal, filfilt_signal, legend="Noisy Signal")
p.line(sine_time, sine_signal, line_width=2, line_color='blue', legend_label='Original signal', line_alpha=0.5)

show(row(p,q,r))

We can see how the application of the one dimensional filter causes a clear misalignment in the signal"s phase. The filtfilt filtering solves this issue, and zeros the phase response to match the original pre-filtered signal.

3 - Application to real life signals

Though the above example has considered a simulated signal, the same procedure can be applied to real life signals , as will be demonstrated in this section

3.1 - ECG signal example

We will use a ECG signal sample, as it is one of the most well known biosignals and has been widely studied before.

In [7]:
# Load of data
data, header = bsnb.load("../../signal_samples/ecg_20_sec_100_Hz.h5", get_header=True)

# The data file is loaded as a Python dictionary, where each key indicates the 
# channel of the acquisition. In the next line we get the channel where the 
# ECG data was acquired.
channel = list(data.keys())[0]

# The header is also in the form of a Python dictionary. The next line saves 
# the sampling frequency of the acquired data
fs = header["sampling rate"]
# Resolution of the sensor (number of available bits)
resolution = header['resolution'][0]

# Signal Samples
signal_raw = data[channel]
time = bsnb.generate_time(signal_raw, fs)

# Let's convert the signal's units, since we know it is a ECG signal
signal = bsnb.raw_to_phy("ECG", "biosignalsplux", signal_raw, resolution, "mV")

After loading the signal and converting it to physical units, we can apply the two types of low-pass filter: normal and filtfilt . Both will be Butterworth filters of order 3 and low-pass frequency of 300 Hz. In this case, the cut-off frequency is not relevant for the analysis of signal, as we only wanted to illustrate the difference of applying the different filters.

In [8]:
# Creating a Butterworth filter with order 3 and low-pass frequency of 300 Hz
lfilter_signal = bsnb.lowpass(signal, 300, order=3)
filtfilt_signal = bsnb.lowpass(signal, 300, order=3, use_filtfilt=True)

The next plots show the results of the application of the different filters. The first shows the original ECG signal without the application of filters. The middle figure shows the original signal and the same signal resulting from the application of a normal low-pass filter. The last figure shows the same, but the signal in dashed line is the result of the application of the filter using filtfilt .

In [9]:
# Let's plot our example
p,q,r = print_filtfilt_plots(time, signal, lfilter_signal, filtfilt_signal)

p.x_range = Range1d(0, 1)
q.x_range = Range1d(0, 1)
r.x_range = Range1d(0, 1)

show(row(p,q,r))

Again, the regular low-pass filter results in a misaligned signal relative to the original signal. However, the filtfilt filter corrects the misalignment and brings the filtered signal back to the original time interval.

3.2 - EEG signal example

The same can be applied to EEG signals, as we will demonstrate in this section.

In [10]:
# Load of data
data, header = bsnb.load("../../signal_samples/signal_sample_single_hub_EEG_2018_7_4.h5", get_header=True)
channel = list(data.keys())[0]

# Sampling frequency of acquired data
fs = header["sampling rate"]
# Resolution (number of available bits)
resolution = header['resolution'][0]

# Signal Samples
signal_raw = data[channel]
time = bsnb.generate_time(signal_raw, fs)

# Let's convert the signal's units, since we know it is a EEG signal
signal = bsnb.raw_to_phy("EEG", "biosignalsplux", signal_raw, resolution, option="V")

This time, we will apply a similar low-pass Butterworth filter of order 3, but using a cut-off filter of 20 Hz.

In [11]:
# Creating a Butterworth filter with order 3 and low-pass frequency of 20 Hz
lfilter_signal = bsnb.lowpass(signal, 20, order=3)
filtfilt_signal = bsnb.lowpass(signal, 20, order=3, use_filtfilt=True)

The next plots show the results of the application of the different filters. The first shows the original EEG signal without the application of filters. The middle figure shows the original signal and the same signal resulting from the application of a normal low-pass filter. The last figure shows the same, but the signal in dashed line is the result of the application of the filter using filtfilt .

In [12]:
# Let's plot our example
p,q,r = print_filtfilt_plots(time, signal, lfilter_signal, filtfilt_signal)

p.x_range = Range1d(7, 8)
q.x_range = Range1d(7, 8)
r.x_range = Range1d(7, 8)

p.yaxis.axis_label = 'Volt (V)'
q.yaxis.axis_label = 'Volt (V)'
r.yaxis.axis_label = 'Volt (V)'

show(row(p,q,r))

The middle figure shows that the regular filter results in a change in the phase of the signal, while the application of the filtfilt corrects the phase of the resulting filtered signal.

3.3 - EMG signal example

In this section, we will show how the same procedure can be applied for EMG signals .

In [13]:
# Load of data
data, header = bsnb.load("../../signal_samples/emg_bursts.h5", get_header=True)
channel = list(data.keys())[0]

# Sampling frequency and acquired data
fs = header["sampling rate"]
# Resolution (number of available bits)
resolution = header['resolution'][0]

# Signal Samples
signal_raw = data[channel]
time = bsnb.generate_time(signal_raw, fs)

# Let's convert the signal's units, since we know it is a EMG signal
signal = bsnb.raw_to_phy("EMG", "biosignalsplux", signal_raw, resolution, "mV")

After loading the signal and converting it to physical units, let"s apply the two types of filter.

In [14]:
# Creating a Butterworth filter with order 3 and low-pass frequency of 150 Hz
lfilter_signal = bsnb.lowpass(signal, 150, order=3)
filtfilt_signal = bsnb.lowpass(signal, 150, order=3, use_filtfilt=True)

The analysis of the plots is analogous to the others: on the left we show the original signal , the middle figure shows the original signal and the result of the application of the normal filter and the right figure shows the result of the application of the filtfilt filter.

In [15]:
# Let's plot our example
p,q,r = print_filtfilt_plots(time, signal, lfilter_signal, filtfilt_signal)

p.x_range = Range1d(1.9, 2.1)
q.x_range = Range1d(1.9, 2.1)
r.x_range = Range1d(1.9, 2.1)

p.y_range = Range1d(-1, 0.39)
q.y_range = Range1d(-1, 0.39)
r.y_range = Range1d(-1, 0.39)

p.legend.location = "bottom_right"
q.legend.location = "bottom_right"
r.legend.location = "bottom_right"

p.yaxis.axis_label = 'milliVolt (mV)'
q.yaxis.axis_label = 'milliVolt (mV)'
r.yaxis.axis_label = 'milliVolt (mV)'

show(row(p,q,r))

Once again, the application of the filtfilt corrects the dephase resulting of the application of the regular filter.

As a conclusion, we have seen that conventional filtering reduces noise in the signal, but has the side effect of delaying the filtered signal"s phase . This poses a problem when we require a high temporal precision for event detection in a time series. Filtfilt handles this problem by applying the same filter forwards and backwards, reducing the time offset to zero .

It is important to mention that filtfilt cannot be used with differentiator and Hilbert FIR filters, or more generally any filter that depends heavily on the signal"s phase.

We hope that you have enjoyed this guide. biosignalsnotebooks is an environment in continuous expansion, so don"t stop your journey and learn more with the remaining Notebooks !

In [16]:
bsnb.css_style_apply()
.................... CSS Style Applied to Jupyter Notebook .........................
Out[16]: